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1  Preface 

This  is  a  written  version  of  a  series  of  invited  lectures  on  differential-algebraic  systems  of 
equations  (DAEs)  at  the  IVth  SERC  Numerical  Analysis  Summer  School  of  Lancaster 
University.  In  line  with  the  aims  of  the  meeting  these  notes  introduce  some  typical 
applications  and  basic  properties  of  DAEs  and  then  present  an  overview  of  recent,  new 
existence  theories  for  such  systems  based  on  differential  geometric  considerations  and 
on  a  numerical  approach  derived  from  these  theories.  In  the  presentation  the  stress  is 
on  general  concepts,  results  and  applications  rather  than  on  detailed  proofs. 

Differential-algebraic  systems  of  equations  (DAEs)  arise  in  many  applications  in 
science  and  engineering.  For  some  examples  we  refer,  for  instance,  to  the  monographs 
13,22]  and  the  many  references  given  there.  Three  typical  applications  are  sketched  in 
Section  2  below.  Over  the  years,  it  has  become  well  known  that  the  solution  behavior 
of  DAEs  may  differ  considerably  from  that  of  standard  ordinary  differential  equations 
(ODEs).  A  valuable  measure  of  the  deviation  of  a  DAE  from  an  ODE  is  the  concept 
of  an  index  which  was  first  introduced  in  [21]  and  has  since  been  formalized  in  various 
ways.  The  index  highlights  also  some  of  the  differences  between  the  existence  behavior 
of  DAEs  and  ODEs  although  it  does  not,  by  itself,  provide  for  any  existence  results. 
In  fact,  up  to  now,  existence  theories  for  nonlinear  DAEs  are  available  only  for  a  few 
selected  classes  of  systems. 

Since  the  solutions  of  any  DAE  are  expected  to  be  smooth  paths  in  some  space 
of  dependent  variables,  we  should  expect  the  equations  to  define  a  dynamical  system 
in  a  suitable  domain  of  that  space.  While  this  connection  with  dynamical  systems 
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is  immediately  obvious  for  ODEs  this  is  certainly  not  the  case  for  DAEs  and  there 
appear  to  be  only  few  studies  that  specifically  address  this  connection  (see  e.g.  37' 
and  38 ).  Some  aspects  of  the  relationship  between  DAEs  and  dynamical  systems  will 
be  discussed  in  Section  3. 

In  a  series  of  papers  38 . 41  '.  and  .36'  a  differential-geometric  approach  has  been 
developed  for  the  analysis  of  the  dynamical  system  underlying  a  DAE  and  for  the 
proof  of  general  existence  and  uniqueness  results  for  such  systems.  In  Sections  4  and  5 
the  results  in  the  two  last -mentioned  papers  are  summarized,  moreover. Section  5  also 
addresses  relations  between  some  of  these  results  and  the  index  concept. 

Comprehensive  introductions  to  numerical  methods  for  DAEs  may  be  found  in  the 
cited  monographs  3.22'.  A  brief  survey  of  some  of  these  methods  is  given  in  Section 
6.  Then  in  that  section  a  new  local  parametrization  approach  for  DAEs  is  presented 
which  derives  naturally  from  the  differential-geometric  existence  theories  and  has  been 
found  to  lead  to  very  promising  methods  for  the  computational  solution  of  higher-index 
DAEs.  For  the  case  of  the  Euler-Lagrange  equations  of  constrained  mechanical  systems 
this  approach  includes  the  socalled  method  of  generalized  coordinate  partitioning  in¬ 
troduced  first  in  :48:. 

2  Model  Problems 

This  Section  provides  three  illustrative  examples  of  practical  applications  leading  to 
differential-algebraic  equations.  As  indicated  before,  there  are  numerous  other  areas 
were  DAEs  occur. 


2.1  Constrained  Dynamical  Systems 

A  major  source  of  DAEs  is  the  kinematic  and  dynamic  analysis  of  mechanical  multi¬ 
body  systems.  This  is  a  venerable  field  of  mechanics  and  we  give  here  only  some  very 
simple  examples  and  refer  for  further  details  to  the  extensive  literature  (see  e.g.  23.49  ). 

Suppose  that,  under  the  influence  of  a  force  Q,  a  particle  with  mass  m  slides  on  a 
two-dimensional  surface  in  R3  specified  by  the  real-valued  equation  $(x)  =  0.  In  order 
for  the  point  to  remain  on  the  surface,  a  constraining  force  must  act  in  the  normal 
direction  of  the  surface.  If  D$(x)  £  ^(R3^1)  denotes  the  derivative  of  $  at  any 
x  €  R3,  then  this  normal  direction  is  given  by  the  vector  D$(x)T  £  R3.  Hence,  by 
Newton's  law  we  obtain  here  the  DAE 

(2.1)  $(x)  =  0,  mx"  -1-  zD$(x)T  =  Q 

where  z  €  R1  specifies  the  size  of  the  constraining  force.  For  example,  suppose  that 
the  surface  is  a  paraboloid  and  that  gravity  is  the  only  force  acting  on  the  mass,  then 
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(2.1)  becomes 

2  2 

”  2! 3 

(2.2)  mi"  —  2cii  =  0 

mij  —  2zx2  —  0 

mij  —  z  =  mg. 

More  generally,  suppose  that  the  vector  i£R"  characterizes  the  configuration  of 
all  bodies  of  a  mechanical  system  and  that  the  kinematic  constraints  acting  on  the 
system  are  modelled  by  the  s-dimensional  (holonomic)  constraint  equations 

(2.3)  $(x,f)  =  0. 

Here  $  is  now  a  mapping  from  Rn  x  R1  into  R*,  1  <  s  <  n  and  t  represents  time.  Then 
the  equations  of  motion  are 

(2.4)  M(x.  t)x"  ~  Dr$(x,  tfz  =  Q{x.  x',t) 

where  M(x,t)  is  the  mass  matrix  and  Q(x,x',t)  the  vector  of  applied  forces. 

As  an  example  consider  a  simple,  planar  "slider  crank”  consisting  of  two  bodies, 
namely,  a  bar  of  length  2  and  a  wheel  of  radius  1  centered  at  the  origin  of  a  (£.77) 
coordinate  system  in  the  plane.  At  one  of  its  ends  the  bar  pivots  around  a  fixed 
point  on  the  circumference  of  the  wheel  while  its  other  end  slides  along  the  £-axis. 
Any  configuration  of  the  system  may  be  characterized  by  the  vector  x  =  (qj.  q2.  £)t 
consisting  of  the  current  coordinate  £  of  the  bar's  sliding  end  and  the  two  angles  Qj  and 
q2  between  the  £-axis  and  the  directions  to  either  end  of  the  bar.  respectively.  Thus, 
if  the  wheel  turns  with  a  constant  torque  r,  then  the  equations  (2.4),  (2.3)  have  here 
the  form 

cos  Qj  +  2  cos  a2  =  £ 
sin  a2  —  2  sin  a2  =  0 

(2.5)  —  z\  sin  ai  +  z2  cos  Qj  =  r 

m£"  —  2 zj  sin  ct2  —  2z2  cos  a2  =  0 

«/2o ij  —  Z\  —  0 

where  Jj  and  J2  are  the  moments  of  inertia  of  the  wheel  and  the  bar,  respectively,  and 
m  is  the  mass  of  the  bar. 

2.2  Electrical  Circuits 

A  second  extensive  source  of  DAEs  is  the  analysis  of  electrical  circuits.  Once  again,  we 
refer  for  details  to  the  literature  (see  e.g.  [10]  )  and  discuss  only  the  basic  ideas  and  a 
simple  example. 
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A  circuit  may  be  considered  as  an  inter-connected  collection  of  electrical  devices, 
such  as.  resistors,  inductors,  capacitors,  sources,  etc.  Its  connection  pattern  is  modelled 

by  a  finite,  directed  graph  Cl  =  (JV.  A)  with  node  set  .V  =  {1.2 . p}  and  branch  set 

A  =  {At . Aq }  C  A  x  A  where  single-node  loops  (z.  z)  —  A  x  A  are  excluded.  Each 

branch  corresponds  to  a  specific  component  of  the  circuit.  As  an  example  consider  a 
graph  with  four  nodes  with  the  following  branches  and  components: 


branch  (1.2) 
branch  (2, 3) 
branch  (1,3) 
branch  (1,4) 
branch  (4. 1) 
branch  (4, 3) 
branch  (3,4) 


linear  resistor  resistance  =  Rx 
voltage  source  voltage  =  u0 
linear  capacitor  capacitance  —  C j 
linear  inductor  inductance  L\ 
linear  capacitor  capacitam  C2 
linear  capacitor  capacitance  =  C3 
linear  resistor  resistance  =  /?2. 


Generally,  the  graph  Cl  can  be  characterized  by  its  (node-arc)  incidence  matrix 
A  E  Rpx<J  with  the  elements 


a*j  = 


+  1 
-1 
0 


if  (z,  k)  E  A  for  some  k  E  N 
if  ( k,i )  E  A  for  some  k  E  N 
otherwise. 


In  our  example  the  graph  underlying  the  circuit  then  has  the  incidence  matrix 


/  1  0  1  1 

-110  0 

0-1-1  0 

\  0  0  0  -1 


-1  0  0  \ 

0  0  0 

0  -1  1 

1  1-1/ 


With  each  branch  A j  =  (z,  k)  of  Cl  two  electrical  quantities  are  associated,  namely 
a  current  y3  and  a  voltage-drop  ur  They  are  connected  by  a  functional  relation 


(2.6) 


=  0,  j  =  1.2 . q, 


the  so-called  branch- characteristic  of  Aj.  The  specific  form  of  (2.6)  depends  on  the  type 
of  the  device  modelled  by  the  branch,  such  as,  for  instance, 


current  source: 
voltage  source: 
ideal  diode: 
linear  resistor: 
voltage  driven  resistor: 
current  driven  resistor: 
linear  inductor: 
linear  capacitor: 


Vj  =  tfr(t) 

Uj  =  ip(t ) 

max(uj,  -y:)  =  0 
uj  =  Ryj 

V:  =  V'(tij) 

=  Myj) 

Cu/  yj 
Lyj  =  u: 
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In  our  example,  the  set  of  branch  characteristics  (2.6)  is  given  by  the  equations 

#i2/i  = 

u2  =  u0 

CiU3'  =  y3 

Ly 4  =  u4 

C2u5'  =  y5 

=  2/6 
#22/7  =  u7. 

Kirchhoff's  first  conservation-law  requires  that  the  (algebraic)  sum  of  the  currents 
on  the  branches  starting  at  a  node  must  equal  the  sum  of  the  currents  on  the  branches 
terminating  at  that  node.  In  terms  of  the  incidence  matrix  .4  this  means  that  a  per¬ 
missible  current  flow  is  characterized  by  any  vector  y  =  (y\.  ■  ■  ■  .yq)T  6  R9  for  which 

(2.7)  Ay  =  0. 

KirchofTs  second  law  specifies  that  the  (algebraic)  sum  of  all  the  voltage  drops 
on  the  branches  of  any  loop  of  fl  has  to  be  zero.  If  we  introduce  the  vector  u  = 

(ui - -  uq)  £  R9  of  voltage  drops,  as  well  as  the  vector  w  =  (wi, ... ,  wp)T  £  Rp  of  all 

nodal  voltage  levels,  then  the  second  law  is  corresponds  to  the  equation 

(2.8)  u  =  ATy  ,  tiq  =  0 

where  the  last  equation  was  introduced  to  fix  the  absolute  values  of  the  nodal  voltages. 

Thus  altogether  (2.6), (2. 7), (2.8)  form  a  DAE  of  2q  4-  p  *  1  equations  in  2q  -  p 
unknown.  The  reason  for  this  difference  is  that  A  does  not  have  full  rank.  If  fi  is  a 
connected  graph  -  which  is  certainly  a  reasonable  assumption  -  then  a  standard  theorem 
of  graph  theory  (see  e.g.  [6])  ensures  that  rankA  =  p  —  1.  Thus  one  of  the  equations 
(2.7),  is  a  linear  combination  of  the  others  and  hence  may  be  dropped. 

The  equations  (2.6),  (2.7),  (2.8)  are  called  a  descriptor  form  of  the  circuit.  There 
axe  many  ways  of  reducing  the  size  of  this  system  but  we  shall  not  enter  into  any  details 
here. 


2.3  Punch- Stretching  of  Sheet  Metal 

We  end  with  a  somewhat  different  example  arising  in  connection  with  sheet  metal 
stamping  processes.  Since  in  this  case  the  formulation  is  somewhat  more  complex,  we 
do  not  include  all  the  details  but  refer  instead  to  the  literature  (see  e.g.,  [8,9]). 

The  processes  to  be  considered  involve  the  deformation  of  a  sheet  of  metal  in  a 
forming  press  with  a  particular  punch  and  die  configuration.  In  order  to  ease  the 
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discussion  we  consider  a  simpler  problem,  namely  the  so-called  hydrostatic  bulge  test 
used  widely  in  metallurgy.  An  initially  fiat  sheet  of  metal  is  clamped  over  one  end  of  a 
cylindrical  chamber  into  which  hydraulic  oil  is  then  pumped.  This  creates  a  hydrostatic 
load  on  the  sheet  and  causes  it  to  bulge  outward. 

In  line  with  the  formulation  presented  in  47'  the  equation  of  virtual  work  for  the 
hydrostatic  bulge  deformation  is  given  by 

(2.9)  h0  [  rSe  dAo  =  f  p6v  dA0. 

J  Ao  J  Ao 

Here  h0,  A0  denote  initial  sheet  thickness  and  surface  area,  respectively.  £  is  the 
radial  distance  to  a  erial  point  on  the  sheet  at  time  zero,  and  v  is  a  volume  measure. 
Moreover,  r  =  ( .  r2  j  is  the  vector  of  the  Kirchhoff  stresses  in  the  radial  and  circumfer¬ 
ential  directions,  respectively,  while  s  =  is  the  vector  of  the  logarithmic  strains 

in  the  corresponding  directions,  defined  by 

-i  =  ln[(l  -  u{)2  -  u£1/2  ,  s2  =  ln'l  -  ^ 

s 

where,  u  and  w  denote  the  radial  and  vertical  displacements  of  the  material  point  whose 
initial  position  is  given  by  (£.0). 

In  addition,  we  have  the  material-dependent  constitutive  equations  in  rate  form. 
These  have  the  generic  form 

(2.10)  t'  =  Ls'  -r  r(r,r),  s' =  g(r,I) 
where 


and  E  and  u  are  Young’s  modulus  and  Poisson’s  ratio,  respectively,  s  is  the  effective 
strain,  and  the  nonlinear  functions  r  and  g  depend  on  the  specific  form  of  the  hardening 
law. 

The  volume  measure  in  (2.9)  is  defined  by 

V  —  1  T  j)(l  +  U€) 

and 

V(t)  =  /  v(0  dAc 

J  Ao 

is  the  volume  of  the  bulge  at  time  t  which  may  be  assumed,  for  instance,  to  satisfy 
V(t)  =  7 1  with  fixed  7. 
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For  the  computation,  we  introduce  finite  element  approximations  of  the  displace¬ 
ments  u.tv,  the  stress  components  rx ,  r2 ,  and  the  effective  stress  i.  Then  the  equation 
(2.9)  of  virtual  work  is  approximated  by  a  nonlinear  equation  of  the  form 

(2.11)  F(x,y.q,t)  =  0 

where  the  vectors  x,y,q  contain  the  approximations  of  (u.ir),  (r1.r2).  and  p.  respec¬ 
tively.  Correspondingly,  the  constitutive  equations  (2.10)  are  approximated  by  a  dif¬ 
ferential  equation  of  the  form 

(2.12)  y'-Lz'  =  My,:) 

-  =  My*z) 

where  the  vector  c  represents  the  approximation  of  the  effective  strain  s.  Thus,  alto¬ 
gether  the  equations  (2.11),  (2.12)  form  a  DAE. 

3  DAEs  and  Dynamical  Processes 
3.1  DAEs  and  ODEs 

The  examples  of  Section  2  provide  an  indication  of  various  possible  forms  of  DAEs.  In 
all  cases  the  differential  equations  and  algebraic  equations  turned  out  to  be  separated. 
Observe  also  that  there  may  be  variables  for  which  no  derivatives  appear  anywhere  in 
the  system.  Of  course,  for  the  computation  various  specific  properties  of  the  equations 
are  of  particular  interest.  For  instance,  it  is  usually  advantageous  when  the  derivatives 
only  occur  linearly,  etc. 

In  some  applications  the  form  of  the  equations  may  vary  in  different  parts  of  the 
space,  and,  in  particular,  there  may  not  exist  a  globally  valid  separation  into  algebraic 
equations  and  differential  equations.  Hence,  in  such  cases,  the  DAEs  have  the  generic 
form  of  an  implicit  differential  equation 

(3.1)  F(x,x',t)  =  0 
which  cannot  be  transformed  into  the  form 

(3.2)  x'  =  f(x,t) 

of  an  explicit  ordinary  differential  equation  (ODE) 

If,  in  (3.1),  F  is  a  sufficiently  smooth  map  from  R2n+1  into  Rn  and  the  derivative 
DpF(x,p,t)  €  T(Rn,Rn)  is  an  isomorphism  at  some  solution  (x,y,t)  of  F(x,y,t)  =  0 
then  the  implicit  function  theorem  guarantees  that  (3.1)  can  be  transformed  locally 
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into  the  form  (3.2).  Hence,  in  our  setting  we  should  assume  that  DpF(x. p.t)  does  not 
have  full  rank.  More  specifically,  we  shall  call  (3.1)  an  implicit  DAE  only  if 

(3.3 )  rank DpF{x,p.t)  =  constant  <  n, 

on  the  domain  under  consideration.  This  constant-rank  assumption  excludes  various 
singular  implicit  equations  (3.1)  with  a  solution  behavior  that  may  differ  radically  from 
that  of  ODEs  or  DAEs  (see  e.g.  ’35]). 

The  existence  and  uniqueness  theory  for  solutions  of  explicit  ODEs  (3.2)  is  a  well- 
developed  subject  (see  e.g.  11]  and  also  the  Appendix).  In  particular,  if  /  :  E  Z 
FT-1  — •  Rn  is  of  class  C1  on  some  open  set  E.  then  we  know  that  for  any  point 
( x0.t0 )  €  E  there  exists  a  C2-solution  x  :  J  C  Rl  —  E  of  (3.2).  defined  on  some  open 
interval  J  containing  t0,  which  satisfies  the  initial  condition  x{to)  —  x0.  Moreover,  any 
two  such  solutions  satisfying  the  same  initial  condition  are  identical  on  the  intersection 
of  their  domain. 

This  local  solvability  result  for  ODEs  does  not  carry  over  directly  to  DAEs.  In  fact, 
consider  the  simple  system 


Xj  =  cos  12 

(3.4)  xj  =  x3 

x'2  =  1, 

and  suppose  that  x  :  J  C  R1  —  R3  is  any  C^-solution  of  (3.4)  on  some  open  interval  J. 
Then,  by  differentiating  X\(t)  —  cos  x2(f)  =  0  with  respect  to  t  and  using  the  differential 
equations,  we  find  that  the  solution  must  satisfy  the  algebraic  condition 

(3.5)  x3  —  sin  x2  =  0, 
whence  necessarily 

(cos  t  \ 
t 

—  sint  / 

Conversely,  (3.6)  does  define  a  C°°  -solution  x  :  R1  — *  R3  of  (3.4).  In  other  words.  (3.6) 
is  the  only  solution  of  (3.4)  and  we  certainly  cannot  prescribe  any  initial  conditions 
other  than,  trivially,  a  point  of  (3.6). 

This  result  indicates  that  for  DAEs  there  may  be  some  ’’hidden”  constraints,  such 
as  (3.5),  which  all  solutions  have  to  satisfy.  As  a  consequence,  there  may  not  exist  ab¬ 
solution  through  every  choice  of  initial  point;  that  is,  only  certain  initial  conditions 
may  be  admissible. 

Obviously,  any  ’’hidden”  constraints  may  be  expected  to  cause  difficulties  during 
the  numerical  solution  of  the  DAE.  This  is  indeed  the  case,  and,  in  fact,  it  is  by 
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now  well  known  that  the  degree  of  difficulty  rises  with  the  number  of  such  additional 
constraints.  This  observation  led  W.Gear  and  L.  Petzold  (  21])  to  introduce  an  index 
which  measures  the  "deviation”  of  a  DAE  from  an  ODE.  We  shall  discuss  later  some 
aspects  of  this  concept;  at  this  moment  it  will  suffice  to  characterize  the  index  of  a 
DAE  loosely  as  the  total  number  of  given  algebraic  and  hidden  constraints  that  are 
needed  to  specify  the  solution  completely.  In  this  sense.  (3.4)  is  a  DAE  with  index  2. 

As  another  example  consider  the  DAE  (2.2)  modelling  the  dynamics  of  a  mass-point 
on  a  paraboloid.  This  system  can  be  written  in  the  first  order  form 

<£(x)  =  x3  —  r’;  —  x3  =  0 

(3.7)  x  =  y 

V  =  ge2  -  zD$(x)T 

where  x,  y  £  R3.  e3  is  the  third  natural  basis  vector  of  R3  and  we  replaced  z:m  by  z. 
By  differentiating  the  algebraic  equation  and  using  the  differential  equations  we  obtain 
as  first  “hidden”  constraint 

(3.8)  2x1y1  -  2x2y2  -  J/3  =  0. 

In  turn,  by  differentiating  (3.8)  we  are  led  to  the  further  constraint 

(3.9)  2{y\  -f  y\)  -  z{\  4-  4(x3  -  x3))  =  0. 

It  is  not  difficult  to  see  that  (3.7)  together  with  these  two  constraints  (3.8).  (3.9) 
completely  specifies  the  solutions.  Thus  in  our  terminology  the  DAE  has  index  3. 

The  two  constraints  have  here  a  simple  geometric  meaning.  In  fact,  (3.8)  requires 
the  velocity  vector  y  to  be  tangential  to  the  paraboloid  while  (3.9)  means  that  the 
constraining  force  zD$(x)T  has  to  balance  the  other  two  forces. 

This  index-result  is  not  restricted  to  the  special  example  (2.2);  in  fact  it  turns  out 
that  all  Euler-Lagrange  systems  (2.3),  (2.4)  have  index  3. 

3.2  Dynamical  Processes 

Ordinary  differential  equations  axe  a  fundamental  tool  in  the  study  of  dynamical  pro¬ 
cesses  that  are  finite-dimensional,  differentiable  and  causal.  By  this  we  mean  processes 
for  which 

(i)  the  states  are  characterized  by  finitely  many  degrees  of  freedom, 

(ii)  the  changes  of  the  states  are  described  by  differentiable  functions,  and 

(iii)  the  future  behavior  is  uniquely  determined  by  the  initial  conditions. 
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Our  examples  suggest  that  differential- algebraic  equations  also  represent  models  of  such 
dynamical  processes. 

The  theory  of  dynamical  processes  has  been  heavily  influenced  by  mechanical  consid¬ 
erations.  Thus,  we  use  a  simple  mechanical  example  to  review  some  of  the  basic  termi¬ 
nology.  Consider  the  motion  of  k  particles  in  R3  and  let  x,  and  y,.  i  =  1.2 . k,  denote 

the  location  and  velocity  of  the  particles  at  time  t.  Then  with  q  =  lx; . x*. )  £  R3;r 

and  p  =  [yi.---.yk)  €  R3k  the  state  of  the  system  is  given  by  (q.p).  The  states  are 
usually  restricted  to  some  specified  subset  5  C  R3fc  x  R3*  -  the  state  space  of  the  system. 

In  many  cases,  the  state  space  is  an  open  subset  of  R3fe  x  R3C  For  instance,  if  no 
two  particles  are  ever  allowed  to  be  in  the  same  place  at  the  same  time,  then 

5  =  €  R3k  x  R3k:  x,  ~  x:  for  i  j.  i.j  =  1. 2, ....  k} 

is  certainly  open. 

On  the  other  hand  -  for  instance  when  there  are  angular  variables  or  constraints 
-  the  state  space  need  not  be  an  open  set.  For  example.  c-\ rider  a  system  of  rigidly 
connected  particles.  Then  the  configuration  vectors  q  have  to  belong  to  the  set 

C  =  {q  6  R3<!;  !j  x,  -  Xj  j!,  =  dji  for  i  ±  j.  i,  j  =  1. 2, . . . ,  k} 

where  dj  are  given  constants.  For  k  >  4  this  is  a  six- dimensional  submanifold  *  of 
R3fc.  This  follows  from  the  well-known  fact  that  the  position  of  a  rigid  body  in  R3  is 
uniquely  characterized  by  the  location  of  one  point  and  the  orientation  of  an  orthonor¬ 
mal  coordinate  system  fixed  within  the  body.  Then  the  state  space  may  be  identified 
with  the  tangent  bundle  C  x  R3k  of  C  and  hence  is  a  12-dimensional  submanifold  of 
R3k  xR* 

In  differential-geometric  terms  the  two  cases  are  not  very  different.  In  fact,  any 
open  subset  of  R3*  x  R3;e  is  a  6A:-dimensional  submanifold  of  that  space  and  hence, 
in  either  case,  the  state  space  is  a  submanifold.  This  agrees  with  the  fundamental 
assumption  introduced  by  H.  Poincare  1880)  that  the  state  space  of  a  mechanical 
system  should  be  t  differentiable  manifold.  Correspondingly,  the  dynamical  system  is 
viewed  as  a  field  of  vectors  on  this  manifold  such  that  a  solution  is  a  smooth  curve 
tangent  at  each  of  its  points  to  the  vector  attached  to  that  point.  We  refer,  e.g.,  to  1 
and  the  historical  references  included  there. 

However,  from  a  computational  viewpoint,  there  is  indeed  a  substantial  difference 
between  the  above  two  cases  which,  in  fact,  reflects  again  the  earlier  indicated  differ¬ 
ences  between  ODEs  and  DAEs.  As  noted  before,  in  the  classical  theory  of  explicit 
ODEs  (3.2)  we  assume  /  to  be  of  class  C1  on  some  open  subset  E  of  R"^1  and,  of 
course,  on  E  the  ODE  induces  the  natural  vector  field 

{x,t)€E  -  ((x,t),/(x,t))  €  TE  =  Ex  Rn+1. 

'Somt  differential-geometric  concepts  and  results  are  collected  in  the  Appendix. 
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Since  the  points  {x.t)  represent  the  states  of  the  system,  this  corresponds  to  the  case 
when  the  state  space  is  an  open  subset.  On  the  other  hand,  as  the  above  example 
of  rigidly  connected  particles  shows,  for  a  DAE  the  state  space  is  expected  to  be  a 
lower  dimensional  manifold.  For  instance,  in  the  trivial  example  (3.4)  this  is  the  one¬ 
dimensional  suomanifold 

{x  €  R3;  —  cos  x2  =  0,  x3  sin  x2  =  0} 

defined  by  the  given  and  hidden  constraints. 

In  the  standard  theory  of  numerical  methods  for  solving  ODEs  of  the  form  (3.2)  it  is 
critical  that  the  domain  E  of  the  right-hand  side  is  open  and  that  there  exists  a  locally 
unique  solution  of  (3.2)  through  each  point  x  of  E.  In  fact,  any  such  method  generates 
a  sequence  {(xk,tk)\k  =  1.2,...}  of  approximating  points  on  the  solution  through  a 
given  initial  point  (io,  to)-  At  best,  we  know  that  these  points  belong  to  some  open 
neighborhood  of  the  exact  solution  contained  in  the  open  set  E  .  Furthermore,  for 
the  step  from  xk  to  x^_i  all  methods  are  designed  to  approximate  the  local  solution 
through  x*.  in  the  sense  that  the  error  between  xk~i  and  this  local  solution  converges 
to  zero  when  the  steplength  tends  to  zero.  It  is  by  no  means  obvious  how  to  extend 
this  approach  to  the  case  when  the  domain  E  is  no  longer  an  open  set  but  some  lower- 
dimensional  submanifold  of  R"  x  R1. 

4  Existence  Theory  for  Implicit  DAEs 

As  noted  earlier,  the  literature  on  DAEs  is  growing  rapidly  but  general  existence  the¬ 
ories  have  only  begun  to  be  developed  relatively  recently.  The  earliest  such  result 
appears  to  be  the  existence  theory  for  gradient  systems 

VxS'(i.y)  =  0,  x  =  f{x,y). 

developed  by  F.Takens  [45]  who  used  the  approximating,  singularly  perturbed  system 
of  differential  equations 


ex'  =  Xxg(x,y)  ,  x  -  f{x,y). 
with  asymptotically  small  t. 

Best  understood  are  probably  the  linear  DAEs  with  constant  coefficients 

(4.1)  Ax'  4-  Bx  =  g(t),  x  €  Rn,  A,  B  €  L( Rn),  rank.4  —  r  <  n 

for  which  existence  results  can  be  proved  by  means  of  the  Kronecker  canonical  form  for 
matrix  pencils  (see,  e.g.,  [16]).  For  a  presentation  of  this  theory  see  e.g.  [21]  or  [22  . 
For  further  references  see  also  [26,46]. 
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In  38  the  indicated  interpretation  of  DAEs  as  dynamical  systems  on  manifolds  was 
used  to  obtain  existence  results  for  semi-explicit  systems 

Fi(x)  =  0,  A(x)x'  =  G(x ). 

These  results  were  generalized  in  41  to  first  and  second  order  systems  of  the  form 

(4.2)  Fa(x)  =  0,  Fa(x.x',z)  =  0. 
and 

(4.3)  F1(x)  =  0,  F2(x,  x\x".  :)  =  0. 
respectively. 

For  general  implicit  equations  (3.1)  -  of  course,  under  the  constant-rank  assumption 

(3.3)  -  local  existence  results  were  first  given  in  ’22],  however,  under  the  restrictive 
condition  that  ker DpF(x.p.t)  is  independent  of  x  and  p.  Finally,  without  such  a 
condition,  a  solution  theory  for  implicit  DAEs  was  presented  in  [36  .  This  theory  will 
be  outlined  in  the  following  sub-section:  for  proofs  and  further  details  we  refer  to  the 
original  article. 


4.1  Local  Theory  for  Implicit  DAEs 

For  ease  of  notation  we  shall  consider  (3.1)  in  the  autonomous  form 

(4.4)  F(x,x')  =  0, 


under  the  three  assumptions 

(Al)  F  :  E  C  Rn  x  Rn  — ♦  Rn  is  C2  on  the  open  set  F; 

(A2)  rankDF(x,p)  =  n,  V  (x,p)  €  F; 

(A3)  rankFpF(x,p)  =  r  <  n,  V  (x,p)  6  E. 

The  condition  (A2)  requires  that  the  equations  (4.4)  are  independent  and  also  implies 
that  F_1(0)  is  an  n-dimensional  C2-submanifold  of  Rn  x  Rn  (see  again  the  Appendix). 
An  instructive  prototype  for  (4.4)  is  the  semi-implicit  DAE 


(4.5) 


F(x,x') 


*!(*)  ^ 

F2(x,x')  ) 


Here  (Al)  holds  if  Fj  :  Ex  —  Rn"r  and  F2  :  Ex  x  Ep  —  Rr,  are  of  class  C2  on 
open  sets  Ex  and  E  —  Ex  -  Fp,  respectively.  Moreover,  (A2)  and  (A3)  are  satisfied  if 
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rankFF^x)  =  n  —  r,  Vx  6  Ex  and  rank DpF2(x,p)  =  r,  t/(x,p)  €  F.  Obviously.  the 
first  of  these  conditions  implies  that 

(4.6)  M  =  {x  €  Ex\  F\{x)  =  0} 

is  an  r-dimensional  C:-submanifold  of  Rn. 

A  (^-solution  of  the  general  equation  (4.4)  is  any  mapping 

x  :  J  -  F",  (x(t).x'(t))  E  E.  F(x(t).x'(t))  =  0,  Vt  6  J 

that  is  of  class  Ck  on  some  open  interval  J  of  R1.  As  noted  in  the  previous  Section,  we 
cannot  expect  that  there  is  a  solution  through  each  point  (x,p)  €  E  and  the  following 
lemma  provides  a  necessary  condition  for  this  to  hold: 

Lemma  4.1  For  (x,p)  6  E  the  conditions 

(4.7)  F(x, p)  =  0,  DxF(x,p)p  €  rgeDpF(x.p). 

are  neccessary  for  the  existence  of  a  Cl  -solution  of  (4-4)  that  passes  through  (x.p). 

For  any  C2-solution  of  (4.4)  this  follows  directly  from  the  fact  that  by  differentiation 
of  F{x{t),  x'(f))  =  0  we  obtain  DxF(x(t).  x'(t))x'{t)  —  DpF(x(t ),  x'(t))x"(t)  =  0  for  all  t 
in  J.  Of  course,  for  C1 -solutions  this  argument  cannot  be  used  and  a  more  subtle  proof 
is  required  (see  (36]). 

For  a  closer  analysis  of  the  set  of  points  characterized  by  the  necessary  conditions 

(4.7)  we  introduce  the  orthogonal  projections 

F,<?  :  E  —  L(Rn,Rn),  P(x,p)Rn  =  rge£>pF(x,p), 

(?(x,p)  =  In  -  P(x.p),  V(x,p)  €  £, 

Because  of  the  constant  rank  condition  (A3)  these  projections  are  C1  -functions  on  E. 
Hence,  also  the  reduced  map 

(4.8)  F  :  E  —  Rn,F(x,p)  =  P(x,p)F(x,p)  4-  Q{x,p)DxF{x,p)p,  (x,p)  €  F, 

is  of  class  C1  on  E.  Then  we  can  show  that  the  set  En  of  all  points  satisfying  the 
necessary  conditions  (4.7)  is  given  by 

(4.9)  En  =  {(x,p)  €  F;  F(x,p)  =  0,F(x,p)  =  0}. 

In  the  special  case  of  (4.5)  the  projections  axe  independent  of  x  and  p.  In  fact  we 
have 
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and  therefore 
(4.10) 


F(x,p) 


DFl(x)p  \ 
F-yix.p)  J 


and 

Ex  =  {(i.p)  €  E{  Fi(x)  =  0,  DFi(x)p  =  0.  F2(x,p)  =  0  }. 
Generally,  the  mapping  (4.8)  defines  the  reduced  equation 


(4.11) 


F(x,  x')  =  0, 


which,  in  essence,  has  the  same  solutions  has  the  original  DAE.  More  specifically  the 
following  result  holds: 


Lemma  4.2  Any  Cl  -solution  of  the  original  equations  (4-4)  solves  the  reduced  equation 
(4-11).  Conversely,  any  C2 -solution  of  (4-11)  that  passes  through  some  point  ofF~l{ 0) 
is  a  C2 -solution  of  (4-4) • 

In  the  special  case  of  (4.5)  when  F2  is  linear  in  p\  that  is,  in  the  case  of  the  DAE 

F(z,x')=(  Affix' -G(x)  )  =0' 

DpF(x,p )  does  not  depend  on  p  and  the  reduced  equation  (4.11)  becomes  the  linear 
equation 

(4.13)  *<*>*'=  (g(i))=0’ 

Suppose  that  the  subset 

(4.14)  M0  =  {x  €  M;B{x)  6  Isom(iZn)} 

of  the  constraint  manifold  (4.6)  is  not  empty.  Then  M0  is  an  r-dimensional  submanifold 
of  Rn  on  which  (4.13)  induces  the  tangential  C^-vectorfield 

x  €  M0  ~  (x,p)  6  TM0,  p  =  S(x)_1  ^  . 

It  is  readily  seen  that  the  integral  curves  of  this  vectorfield  are  exactly  the  solutions  of 
(4.13)  and  therefore,  by  Lemma  4.2,  also  of  the  DAE  (4.12).  This  corresponds  to  the 
approach  used  in  [38]  to  develop  an  existence  theory  for  DAEs  of  the  form  (4.12). 

For  the  general  DAE  (4.4),  we  proceed  analogously  and  assume  that  the  set 

EA  =  {(x,p)  €  En\  DpF(x,p)  €  Isom(Rn)} 


(4.15) 
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is  non  empty.  Clearly,  by  continuity  EA  is  (relatively)  open  in  £\-.  Moreover,  by  the 
implicit  function  theorem  it  follows  that  locally  in  some  open  neighborhood  of  any 
point  (xo.po)  £  Ej i  the  reduced  equation  (4.11)  can  be  transformed  into  an  explicit 
ODE.  Thus  by  applying  the  standard  ODE-theory  we  can  now  prove  the  following 
local  existence  result  for  the  original  DAE  (4.4): 

Theorem  4.1  Given  t0  £  Rn.  consider  the  initial  value  problem 

(4.16)  F{x.x')  =  0.  x(t0)  =  x o,  x'(t0)  =  po 

•mder  the  assumptions  (Al.2.3).  If  a  Cl-solution  of  (4-16)  exists,  then  {x0.p0)  £  E\. 
Conversely,  for  any  (loiPo)  £  EA  there  exists  a  C1  -solution  of  (4-16)  which  is  unique 
on  some  sufficiently  small  interval  J  containing  t0.  Moreover,  this  solution  is  actually 
of  class  C~  on  J . 

The  set  EA  of  (4.15)  has  a  manifold  structure.  This  is  self-evident  in  the  case  of 
(4.5)  where,  obviously,  the  derivative  of  the  mapping 

(x’rt iE  ~  {  FCr)  )  e R*~ 

has  full  rank  for  any  point  of  EA  and  hence,  EA  is  either  empty  or  an  r-dimensional 
C1 -submanifold  of  R"  x  Rn.  But  the  result  also  holds  in  general: 

Lemma  4.3  The  set  EA  C  E  of  admissible  initial  points  of  (4-4)  w  either  empty  or 
an  r -dimensional  Cl  -submanifold  o/Rn  x  Rn. 

From  the  implicit  function  theorem  it  follows  that  for  any  point  (x0.p0)  £  EA  there 
exist  an  open  neighborhood  U  =  Sx  x  Sy  C  E  and  a  unique  Catnapping  V  :  Sx  —  Sy 
with  r]{xo)  =  po  such  that  F\x,  p)  =  0  for  (x,p)  £  U  if  and  only  if  p  =  jj(x).  Since  EA 
is  open  in  £y  we  may  assume  that  U0  =  U  n  EA  =  U  H  E^-.  Let  II  :  Rn  x  R"  —  Rn  be 
the  projection  onto  the  first  factor.  Then,  the  result  means  that  the  restriction  II  T'0  is 
a  CCdiffeomorphism  from  U$  onto  nLr0.  Hence 

M  =  SX  n  IK70 

is  an  r-dimensional  submanifold  of  Rn  and 

x  €  M  *— ►  (x,p)  €  TM,  p  =  p(x), 

is  a  tangential  C'-vectorfield  on  M  for  which  it  can  be  shown  that  the  integral  curves 
axe  exactly  the  solutions  of  (4.4)  in  Eq.  In  other  words,  the  following  result  holds: 
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Theorem  4.2  With  the  above  terminology,  x  :  J  —  Rn  is  a  C1  -solution  of  F(x.  x')  =  0 
satisfying  |x(£),x'(t))  G  U  fort  G  J  if  and  only  if  x(t)  G  M.  x  '(t)  =  T](x(t)).  G  J- 

Thus,  as  expected,  the  DAE  (4.4)  is  locally  equivalent  to  an  explicit  ODE  on  an 
r- dimensional  submanifold  of  Rn. 

The  results  in  this  section  can  be  extended  easily  to  the  general  nonautonomous 
case 

F(t,  x,  x')  =  0. 

In  fact,  as  usual  we  can  transform  this  problem  into  autonomous  form  by  :ntroducing 
the  mapping 


G  ■■  Rnil  x  R-1  -  Rn~\  G((t,  x).  (r.p))  =  FT[t  xlp) 

Then  under  the  required  smoothness  assumptions  the  necessary  conditions  (4.7)  for  G 
assume  the  form 

F(t,x,p)  =  0,  DtF{t,x,p)  ~  DxF{t.x,p)p  G  rgeDpF(t,x.p). 

Similary,  we  can  derive  the  form  of  the  reduced  mapping  and  of  the  set  EA  of  admissible 
initial  points  (see  [36]). 

4.2  Globalizations 

Consider  again  the  implicit  DAE  (4.4)  under  the  assumptions  (A  1.2.3).  As  we  saw  in 
Theorem  4.1,  for  any  (xo,po)  in  the  set  EA  of  (4.15)  the  initial  value  problem  (4.16) 
has  a  local  (^-solution  on  a  sufficiently  small  open  interval  J  —  (a,  b )  G  Rn  containing 
t0-  As  in  the  standard  ODE-theory,  under  appropriate  conditions  these  local  solutions 
can  be  continued. 

In  [36]  the  following  basic  continuation  result  was  proved: 

Theorem  4.3  Suppose  that  EA  =  En- 

( i )  If  for  some  e  G  (0,  b  —  a)  the  set  {p  G  Rn;  p  =  x'(t),  b  —  e  <  t  <  b}  is  bounded  then 
lima._b_  x(t)  =  X{,  exists. 

(n)  //limt_b_  z(£)  =  Xb  exists  and  for  some  sequence  {tj,}  G  J  with  limk—actk  =  b, 
the  sequence  {x'(tfc)}  has  an  accumulation  point  p *  for  which  (xj >,p*)  G  E  then 
limt_i>_  x'(t)  =  p*  and  hence,  for  b  <  oo  the  solution  can  be  continued  to  the 
right. 
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An  analogous  result  holds  for  the  left  endpoint.  Thus  any  local  solution  of  (4.16) 
can  be  extended  to  a  maximal  open  interval  Jm  =  (a*.fe*),  —oc<a’<b'<oc. 

For  an  explicit  initial  value  problem 

x  =  f(x).  x{0)  =  xq,  x  £  Rn 

with  smooth  /  on  all  of  Rn.  any  maximally  extended  solution  for  which  lim;_9_  x<  t ) 
exists  has  bounded  derivatives  x'(t)  near  b.  This  does  not  carry  over  to  DAEs  as  the 
following  result  shows  (see  again  36  ): 

Theorem  4.4  Suppose  that  E  =  Rn  and  EA  =  E y.  Ifbm  <  oe,  then  x'(f )  is  unbounded 
on  the  interval  (bm  —  e,bm )  £  J *  for  all  sufficiently  small  e  >  0.  Hence.  z/limt_(,_  x{t J  = 
xb  exists,  then  lim^b-  |jx'(t)j!  =  oc. 

Again,  an  analogous  result  holds  for  the  left  endpoint. 

Theorem  4.2  showed  that  the  implicit  DAE  (4.4)  is  locally  equivalent  to  a  vector 
field  on  some  r- dimensional  submanifold  of  Rn.  These  local  vectorfields  can  be  extended 
by  applying  the  theory  of  covering  spaces.  This  was  first  used  in  )41  in  connection  with 
the  DAEs  (4.2)  and  (4.3)  and  then  extended  to  the  implicit  case  in  [36\ 

We  sketch  only  briefly  the  general  approach.  Clearly,  the  local  result  shows  that 
the  restriction  II|£4  is  a  local  homeomorphism  between  £4  and  II£'J4.  Let  E\  be  some 
non-empty,  arc-connected  subset  of  £4  for  which  (£^,II^),  with  11^  =  II  E\.  is  a 
covering  space  of  II£^.  In  other  words,  each  point  z  £  II £^  is  assumed  to  have  an 
open,  arc-connected  neighborhood  U  such  that  each  arc-component  of  (11^ )-1  U  is  not 
empty  and  is  mapped  topologically  onto  U  by  II£^.  Often  EmA  =  £4  can  be  used  here. 
This  is  certainly  the  case  when  for  fixed  x  £  II£4  there  are  only  finitely  many  p  with 
(z,  p)  £  £4.  For  instance,  this  holds  for  the  semi-linear  DAE  (4.5).  In  general,  it  is 
always  possible  to  choose  £4  as  the  closure  of  a  non-empty,  pre-compact,  (relatively) 
open,  and  arc-connected  submanifold  of  £4. 

For  any  given  (io,Po)  6  £4  let  now  Af*  be  a  non-empty,  (relatively)  open,  simply 
connected  subset  of  II£^  that  contains  x0.  For  any  z  £  A/*  choose  a  path  £  :  J  —  A/* 
which  connects  xo  with  z.  Then  there  exists  a  unique  lifting  £*  :  J  —  £4  with  initial 
point  (xoiPo)  for  which  II^*  =  £.  This  lifted  path  has  a  unique  endpoint  (x.p)  in 
£4  because  all  paths  in  M *  between  x0  and  x  are  homotopic.  Since  x  was  arbitrary 
in  A/*  our  local  result  can  now  be  used  to  prove  that  Af*  indeed  is  an  r-dimensional 
submanifold  of  Rn  and  that  the  DAE  (4.4)  induces  a  tangential  vector  field  on  A/*  for 
which  all  integral  curves  in  Af*  are  solutions  of  (4.4). 
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5  DAEs  with  Higher  Index 

5.1  Linear,  Constant  Coefficient  DAEs 

As  mentioned  before,  the  linear  problems  with  constant  coefficients  (4.1)  probably 
represent  the  most  extensively  studied  DAEs  in  the  literature.  For  sufficiently  smooth 
g  the  necessary  condition  (4.7)  turns  out  to  have  the  form  Bx  —  g'{t)  -  rge.4  and  it  is 
easily  seen  that  E 4  =£  0  exactly  if 

(5.1)  Au  =  0  and  Bu  €  rgeA  imply  u  =  0. 

The  cited  existence  theory  for  (4.1)  ensures  solvability  if  and  only  if  the  matrix 
pencil  [A.  B)  is  regular;  that  is,  if  there  is  some  A  6  R1  such  that  B  —  A.4  is  invertible. 
A  central  concept  in  the  solvability  theory  of  (4.1)  is  its  index  which  is  defined  to  be 
the  index  of  the  coefficient  pencil  assumed  to  be  regular.  For  any  regular  pencil  (A.B) 
let  A  be  such  that  B  —  A  A  6  Isom(R”),  then  the  index  of  the  pencil  is  the  smallest 
integer  k  such  that 

ker[(f?  —  AA)"1^  +  =  ker[(B  —  A.4)_1.4j 

It  can  be  shown  that  k  is  finite  and  independent  of  the  choice  of  A.  (see  e.g.  22  ). 
and  it  is  also  readily  seen  that  k  =  0  if  and  only  of  A  is  invertible:  that  is.  if  (4.1)  is 
equivalent  with  an  explicit  ODE. 

In  order  to  relate  the  theory  of  Section  4  to  this  index- concept,  let  P  £  I(Rn)  be 
the  orthogonal  projection  onto  rge.4  and  set  Q  =  In  —  P.  As  before  we  differentiate  the 
DAE  (4.1)  and  then  multiply  the  resulting  equation  Ax"  —  Bx'  =  g’{t)  by  Q  in  order  *o 
remove  again  the  second  derivative  of  x.  Together  with  the  projection  of  the  original 
equation  onto  rgeA  this  produces  the  reduced  DAE 

(5.2)  .4jx'  -  Blx  =  gi(t),  Ax  =  P.4  QB,  Bl  =  PB.  gx(t)  =  Pg{t)  -  Qg'[t). 

Even  without  recourse  to  the  earlier  theory,  it  is  readily  checked  that  (5.1)  is  equiv¬ 
alent  with  Ai  €  Isom(Rn)  and,  hence,  that  when  (5.1)  holds  then  (5.2)  can  be  trans¬ 
formed  into  an  explicit  ODE. 

Suppose  therefore  that  A\  is  singular.  Then  we  may  apply  the  same  procedure 
repeatedly,  as  often  as  necessary,  to  obtain  a  sequence  of  DAEs  of  the  form 

(5.3)  '  AjX1  4-  B:x  =  gj{t),  .7  =  0,1,..., 

where  A:,  B:,  g:  are  specified  recursively  by  Aq  =  .4,  B0  =  B,  g0  =  /  and 

(5.4)  Aj+i  =  P :Aj  +  QjBj ,  Bj+ 1  =  PjBj , 

9j+ l(f)  =  Pj9j{t)+  Qj9j(t),  ;=0,1 . 
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while  P:  denote  the  orthogonal  projections  onto  rge.4j  and  Q:  =  In  —  Py  The  process 
stops  with  the  smallest  integer  k  such  that  .4*  is  invertible. 

The  following  result,  proved  in  36],  shows  that  this  integer  k  is  exactly  the  index 
of  the  DAE: 

Theorem  5.1  If  the  matrix  pencil  ( A.  B )  is  regular  and  rank.4  <  n  (so  that  k  >  l) 
then  k  =  k  <  oc.  Conversely,  if  k  <  oc  then  ( A.  B )  is  regular  and  k  =  «. 

In  36  it  was  also  shown  that  the  theory  of  Section  4  provides  all  the  solutions  of 
(4.1)  provided  only  that  g  is  smooth  enough. 

5.2  Nonlinear  Problems  with  Higher  Index 

The  discussion  of  the  previous  section  suggests  that  we  may  proceed  analogously  when 
the  set  Ea  is  empty  for  the  general  implicit  initial  value  problem 

(5.5)  F(x.x')  =  0.  x(0)  =  x0.  x'(0)  =  p0. 

The  first  step  in  the  construction  of  a  sequence  of  problems  corresponding  to  (5.3).  (5.4) 
was  already  done  in  Section  4.  In  fact,  we  differentiated  the  DAE  and  then  applied  the 
projections  P  and  Q  to  obtain  the  reduced  equation  (4.11). 

Our  sufficient  condition  is  that  DpF(x,p)  is  invertible  at  the  given  initial  point 
(i0.po)  €  E\.  If  this  sufficient  condition  does  not  hold.  then,  as  in  the  linear  case,  it 
is  natural  to  construct  recursively  the  sequence  of  mappings 

(5.6)  F°  =  F,  F1  =  F, 

FJ+l  =  Pj{x,p)FJ(x,p)  t  Qj{x,p)DxF3(x,p),  j  =  0, 1,. . . 

where  P:  again  is  the  orthogonal  projection  onto  £>PFJ  and  Q:  =  /„  —  Pr  The  process 
is  repeated  until  DpFk  is  invertible  at  the  point  under  consideration. 

As  before,  one  might  consider  calling  this  integer  k  the  local  index  of  the  problem 
at  the  particular  point.  However,  the  situation  differs  here  in  a  critical  aspect  from 
that  of  the  linear  case.  The  theory  of  Section  4  can  be  applied  to  the  map  Fk  only 
if  the  conditions  (A  1,2, 3)  are  valid  for  all  the  maps  (5.6)  in  some  neighborhood  E0  of 
the  point  (i0,Po)-  In  particular,  we  require  rankDpFJ(i,p)  to  be  constant  in  such  a 
neighborhood  for  the  projections  P:,  Qj  and  hence  F:  to  be  of  class  C1.  In  addition, 
we  also  need  the  three  conditions  to  conclude  from  the  non- singularity  of  DpFk  that 
the  system  Ffc(x,x')  =  0  can  be  tranformed  locally  into  an  explicit  ODE  and  that  the 
original  problem  (5.5)  has  a  unique  solution.  As  mentioned  earlier,  the  existence  theory 
for  (5.5)  changes  considerably  when  the  constant-rank  condition  is  violated  (see  35]). 

In  line  with  this,  the  problem  (5.5)  will  be  defined  to  have  local  index  k  at  (x0,Po) 
if  there  is  some  open  neighborhood  Eo  of  that  point  such  that  for  j  =  0, 1 . k  —  1 
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the  mappings  (5.6)  satisfy  the  conditions  (Al.2.3)  and  DpFk(x0,p0)  is  invertible.  Note 
that  this  index  does  not  merely  depend  on  information  at  the  given  point  and  hence 
has  a  global  nature.  Obviously,  the  theory  developed  in  Section  4  assumes  that  the 
implicit  problem  (4.4)  has  index  one. 


5.3  Semi-Implicit  Problems  with  Higher  Index 

The  recursive  analysis  outlined  in  the  previous  section  is  a  powerful  tool  for  the  theo¬ 
retical  study  of  higher  index  problems.  But  for  specific  classes  of  equations  it  is  often 
easier  to  derive  existence  and  uniqueness  results  directly.  As  an  example  of  this  we 
consider  in  this  section  the  special  problems  (4.2)  and  (4.3)  both  of  which  have,  in 
general,  index  higher  than  one. 

More  specifically,  for  the  system 


(5.7) 


F(x.x',z ) 


Fi(x)  \ 

F2(x.x\z)  ) 


suppose  that 

(Bl)  Fl  :  Ex  €  Rn  -  R*  is  C2, 


(B2)  F2  :  E2  =  Ex  x  Ep  x  Ez  €  Rn  x  Rn  x  Rm  —  Rr  is  C\ 


(B3) 


DFl(x) 


0 


DpF2(x,p.z)  DzF2(x,p.  z) 


)s 


Isom(Rm~n),  V(x,p.z)£E2 


where  s<n<s-rr  =  n~-  rr  and  Ex.  Ep  C  Rn.  Ez  C  Rm  are  non-empty,  open  sets. 

Evidently,  these  conditic. i  imply  that  (Al.2,3)  are  satisfied  for  (5.7).  The  assump¬ 
tion  (B3)  is  often  called  the  index-two  condition  since  the  index  of  (5.7)  can  be  at  most 
two.  Of  course,  in  the  degenerate  case  s  =  m  =  0,  (B3)  implies  that  (5.7)  can  be 
transformed  into  an  explicit  ODE  and  thus  has  index  zero.  Moreover,  it  is  easily  seen 
that  for  s  =  0,  m  ^  0  or  m  =  0,  the  index  is  one. 

From  (B3)  it  follows  that  rankDFj(x)  =  s  whence  each  member  of  the  family  of 
sets 

(5.8)  Mb  =  {x  €  Ex\  Fi(x)  =  6},  fee  Fi(E«) 


is  an  (n  —  s)-dimensional  C2-submanifold  of  Rn.  Evidently,  any  x0  €  Ex  belongs  to  the 
unique  constraint  manifold  (5.8)  specified  by  fe  =  Fj ( x0 ) .  Moreover,  for  a  solution  of 
(5.7)  to  pass  through  this  point  we  must  have  x0  €  M0. 

Let  x  :  J  €  Ex  be  any  C^-path  defined  on  some  open  interval  J  6  R1.  If  x  is  a  p 
on  that  is,  if  Fi(x(t))  =  0,  Vt  e  J ,  then  necessarily 


(5.9) 


DFi(x(t))x'{t)  =  0, 
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that  is.  t  £  J  —  (x(t),x'(t))  has  to  be  a  path  on  the  tangent  bundle  of  A/0.  The 
explicit  use  of  the  differentiated  constraint  equation  DFx(x )x'  =  0  is  a  basic  step  in 
the  so-called  index-reduction  technique  for  rewriting  the  DAE  as  a  lower  index  system 
(see  18  ).  But  (5.9)  can  also  be  interpreted  in  another  way.  In  fact,  if  (5.9)  holds  for 
some  Cx-path  x  :  J  —  Ex.  then  it  follows  from  the  integral  mean  value  theorem  that 
Fi(x(t))  =  Fj(x(to))  for  any  fixed  t0  £  J.  Hence  t  £  J  —  (x(t).x'(f))  is  a  path  on 
the  tangent  bundle  of  the  manifold  Mb  specified  by  b  =  F1( x0).  This  suggests  that  we 
imbedd  (5.7)  into  the  family  of  DAEs 


(5.10) 


Fj(x)  =  b,  b  £  FX(EX) 
F2(x.x',z)  =  0 


indexed  by  the  vectors  of  Fi(Ex). 

Let  x  :  J  — *  Ex ,  z  :  J  —  E.  be  a  Cx-solution  of  a  member  of  (5.10).  Then  for  any 
point  (x0.po--o)  =  (x(t0),  x'(t0).  z(t0)),  t0  £  J  the  value  b  =  Fx( x0)  uniquely  specifies 
the  particular  DAE.  But,  in  addition,  (x0,p0,-o)  must  satisfy  DFi(x0)p0  =  0.  as  well 
as  F2(xo,Po--o)  =  0.  This  suggests  the  definition  of  the  Cx-map 


(5.11) 


H  :  E2  —>  R*  x  Rr.H{x,p,z)  = 


DFj(x)p 

F2{x,p}z) 


V(x,p,;)  6  E2 


the  initial  data  map  of  the  family  (5.10). 

For  any  given  (x,p,z)  £  E2  the  derivative  DP  ,H  of  H  with  respect  to  (p.  z)  is 
exactly  the  linear  operator  in  condition  (B3).  The  nonsingularity  of  Dp  zH  implies  the 
following  result: 


Lemma  5.1  For  any  ( xo,po ,  2o)  €  K  we  have  an  open  neighborhood  V  =  Sx  x  Sp  x  Sz 
in  E2,  and  unique  Cl-maps  V  :  Sx  -+  Sp,  C  ’■  Sx  ->  Sz,  with  tj(x0)  =  po,  C(*o)  =  m)- 
such  that  for  any  x  £  Sx  the  only  solution  (p,  z)  £  Sp  x  Sz  of  H(x,p,z )  =  0  is 

p  =  p(x),  2  =  C(*)- 


Thus,  on  the  open  neighborhood  Sx  C  Ex  of  x0, 

(5.12)  7T  :  Sx  —  TSX,  n(x)  =  (x,p(x)),  x  £  Sx. 

constitutes  a  CLvectorfield.  Since  DpzH  £  Isom(Rn+m),  the  mapping  H  is  a  submer¬ 
sion  and  hence  the  solution  set 


(5-13)  K  =  {(x,p,z)  £  E2;  H(x,p,z)  =  0} 

is  an  n-dimensional  CTsubmanifold  of  R2n+m.  This  manifold  turns  out  to  be  state 
space  of  the  family  of  DAEs  (5.10).  In  fact,  it  can  be  shown  that  7r  is  tangential  to  the 
constraint  manifold  through  x;  that  is, 

?r(x)  £  TcAfy,  b  =  Fi(x)  Vx  £  Sx, 
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and  that  for  any  solution  x  :  J  —  Sx  of  the  explicit  ODE 
(5.14)  x  =  77(1). 

we  obtain  the  C^-solution  t  £  J  —  (x(t), £(x(t)))  of  the  member  of  (5.10)  specified  by 
b  =  Fi(x[t0))  for  arbitrary  fixed  t0  €  J- 

Thus,  from  the  standard  existence  theory  of  initial  value  problems  for  ODEs  we 
obtain  the  following  result: 

Theorem  5.2  Suppose  that  the  conditions  (Bl.2.3)  hold  and  that  K  is  non-empty. 
Then  any  point  (xo.poi  :o)  €  K  has  an  open  neighborhood  V  =  Sx  x  Sp  x  Sz  Z  E2  such 
that  for  any  xc  G  Sz  there  • s  exactly  one  point  (xc,pc.zc)  G  A’  ^  U.  Moreover,  for  any 
xc  €  Sx  there  exists  a  unique,  maximally  extended  C1  -solution  x  :  J  —  Sx.  z  :  J  —  5., 
on  some  open  interval  J  with  0  €  J,  of  the  DAE  (5.10)  specified  by  b  =  Fi{xc)  which 
satisfies  x'(J)  C  Sp  and  the  initial  conditions  x(0)  =  xc,  x'(0)  =  pc.  r(0)  =  rc. 

We  refer  to  [41]  for  a  globalization  of  this  result  based  on  the  techniques  from  the 
theory  of  covering  spaces  mentioned  at  the  end  of  Section  4. 

The  result  extends  to  the  second  order  DAEs  (4.3).  In  analogy  with  (Bl.2.3)  we 
suppose  that  the  problem 


(5.15) 


F(x,x',x",z)  =  ^  f 


Fi(x) 

2{X,X'.X".Z) 


satisfies  the  conditions: 


(Cl)  F1:Exe  Rn  -  R*  is  C3, 

(C2)  F2  :  E2  =  Ex  x  Ey  x  Eq  x  Ez  6  R3n^m  —  Rr  is  C\ 


DFi(x)  0 

DpF2(x,y,q,z)  DxF2{x,y,q,z) 


G  Isom(Rm*n)  for  each  (x,p,  q.  x)  G  E2. 


where  s<n<s  +  r  =  n  +  m  and  Ex,  Ev  Eq  C  Rn,  Ez  C  Rm  again  are  non-empty, 
open  sets. 

It  is  natural  to  reduce  (5.15)  to  a  first  order  system  by  introducing  a  new  variable  y 
and  adding  the  equation  x'  =  y.  Then  it  turns  out  that,  with  (x,y)  as  new  differential 
variable,  the  resulting  system  constitutes  a  DAE  of  the  form  (5.7)  for  which  (B1.  and 
(B2)  tire  valid.  However,  in  general,  (B3)  does  not  hold  which  is  hardly  surprisin..  since 
we  should  expect  (5.15)  to  induce  local  second  order  vectorfields  instead  of  the  local 
first  order  fields  (5.12). 


5  DAES  WITH  HIGHER  INDEX 


23 


If  x  :  J  —  Ex  is  a  C2-path  on  Mb  for  some  b  G  Fi{Ex):  that  is.  if  F1(x(t))  =  b  for 
t  G  J.  then  for  all  t  G  J  the  we  must  have 

DFi(x(t))x'(t)  =  0, 


as  well  as 

(5.16)  D2Fi(x(t))x"(t)  —  DFi(x(t))(x'(t).x'(t))  =  0. 

This  shows  that  t  G  J  ~  ((x(t),  x'(f)),  (x'(t),  x"(f)))  is  a  path  on  the  second  tangent 
bundle  T2Mb  of  Mb-  Conversely,  bv  the  integral  mean-value  theorem  we  obtain  the 
following  result: 

Lemma  5.2  Let  x  :  J  —  Ex  be  any  C2-path  that  satisfies  (5.16).  If  there  exists  a  to  in 
J  such  that  DFi(x(t0))x'(t0)  =  0  and  therefore  (x(t0),x'{t0))  G  TMb  forb  =  Fi(x(to)). 
then  ((x(t)-x'(t)),  ( x'(t),x"(t )))  €  T2Mb.  for  all  t  6  J. 

As  in  the  first  order  case  this  suggests  that  we  imbedd  (5.15)  into  the  family  of 
DAEs 


(5.17) 


Fi(x)  =  6,  beF,(Ex) 
F2(x,x\x",z)  =  0 


and  that  we  define  the  CMnitial-data  map: 


(5.18)  H:E2 

H{x,y,q,z) 


R*  x  Rr 


(  DF1(x)q  +  D2Fl(x)(y,y) 
V  F2  (x,y,q,z) 


,V(x,y,q,z)  G  E2. 


By  (C3)  we  have  Dq  zH(x,y,q,  z)  €  Isom(R*+r.  Rn+m)  for  (x,y,q.z)  G  E2,  and  hence 
H  is  a  submersion  and  the  solution  set 


(5.19)  K  =  {(x,y,  q,  z)  G  E2\  H{x,y,q,z)  =  0} 

is  a  2n-dimensional  C1  -submanifold  of  R3n+m  and  the  following  result  holds: 

Lemma  5.3  For  any  (xo,  yo,  <?o>  2o)  €  K  there  exists  an  open  neighborhood  U  =  5X  x 
Sv  x  S,  x  Sz  in  E2,  and  unique  Cl-maps  77  :  50  =  5E  x  Sy  — ►  Sp,  C  :  Sx  —  5J;  twth 
tj(x o,yo)  —  qo,  C(xo,yo)  —  such  that  for  any  given  ( x,y )  G  So  the  only  solution 
{q,z)  G  5P  x  Sz  of  H(x,  y,  q,  z)  =  0  is  given  by  q  =  rf(x,y),  z  =  C(*,y). 

Using  this  lemma  we  can  now  define  on  the  open  neighborhood  So  of  (xo,jio)  the 
second-order  CUvectorfield 


(5.20)  ir  :U  <=TSX  ^T2Sx,  ir(x,y)  =  ((x,y),  (y,  77(1,  y))),  V(x,y)  G  S0. 
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Then,  for  any  (x.y)  t  50  such  that  DFfix)y  =  0.  it  follows  that  ir(x.y)  €  T'Mi,  for 
b  =  F-jx).  Moreover,  if  (x,y)  :  J  —  S0  is  any  solution  of  the  explicit  ODE-system 

(5.21)  x  -  y,  y  =  fix.  y), 

satisfying  DFi{x(t0))y(t0)  =  0  for  some  t0  £  d,  then  t  —  (x(t).  i,‘(x(t)) )  is  a  C'-solution 
of  (5.17)  for  b  =  F1(x(t0)). 

Thus  the  standard  solution  theory  provides  here  the  following  local  existence  result: 

Theorem  5.3  Suppose  that  the  conditions  (Cl. 2. 3)  hold  and  that  K  .  •  non-empty. 
Then  any  (x0.  yo- Poi  Ai)  £  F  has  an  open  neighborhood  U  —  *  x  5V  •  Sq  x  5.  m  £2 
such  that  for  any  (xc,yc)  €  S0  =  SxxSy  there  is  exactly  one  po .  -  (xc.yc.pc.  zc)  c  K~~i  . 
Moreover,  for  any  (xc,yc)  6  So  with  DFi(xc)yc  —  0  there  ex; it s  a  unique,  maximally 
extended  C1  solution  x  :  J  £  Sx,  z  :  J  £  Sz  of  (5.17)  for  b  =  F\(xc)  on  some  open 
interval  J  C  Rl  containing  the  origin  which  satisfies  x'(J)  C  Sy.  x"(J)  C  Sq  and  the 
initial  condition  x(0)  =  xc,  x'(0)  =  yc .  x"(0)  =  qc,  r(0)  =  zc. 

Once  again  covering-space  theory  can  be  used  to  globalize  this  result  (see  41  J. 
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6.1  Application  of  ODE  Methods 

The  most  frequently  used  approach  to  the  computational  solution  of  DAEs  is  the  ap¬ 
plication  of  standard  ODE  methods.  This  idea  appears  to  be  due  to  W.Gear  .17:  who 
proposed  the  use  of  backward- difference  (BDF)  methods  in  the  form  developed  for  stiff 
ODEs. 

Briefly,  in  an  m-step  BDF- method  the  derivative  x'  of  the  unknown  function  at 
the  time  tk,  k  >  m,  is  approximated  by  the  derivative  of  the  interpolation-polynomial 
through  (xk,tk)  and  m  earlier  computed  points  i  =  l,2,...,m.  Hence,  in 

the  case  of  the  implicit  initial  value  problem 

(6.1)  F(t,x,x')  =  0,  x(t0)  =  x0,  x'(to)  =  po 

the  determination  of  x requires  the  solution  of  the  nonlinear  system  of  equations 

1  m 

(6.2)  F(tki  Xfc,  &k,j xk— t )  =  Q- 

t=0 

Here  a^i,  i  =  0, 1, . . . ,  m  are  the  coefficients  of  the  BDF  formula  at  the  Ar-th  step  which, 
of  course,  depend  on  k  unless  the  stepsizes  hi  =  U  —  U. x  remain  constant.  For  m  <  7 
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the  m-step  BDF  methods  are  known  to  be  stable  when  applied  to  ODEs  and  hence  we 
assume  from  now  on  that  1  <  m  <  6. 

When  (6.1  I  represents  a  D  ^E:  that  is.  when  the  constant-rank  condition  (3.2  i  holds, 
then  the  validity  and  performance  of  the  process  depends  on  several  factors.  In  par¬ 
ticular.  the  initial  value  problem  (6.1)  has  to  possess  a  solution  and  for  m  >  1  the 
required  m  additional  starting  points  have  to  approximate  this  solution.  Moreover,  at 
each  step  the  nonlinear  system  (6.2)  has  to  have  a  feasible  solution  which  is  computable 
by  a  suitable  iterative  process  such  as  some  form  of  Newton's  method.  In  general,  the 
answers  to  these  questions  depend  strongly  on  the  index  of  the  DAE. 

For  simplicity,  we  restrict  ourselves  here  to  the  semi-implicit  equation  (4.5):  tha*  is. 
to  the  initial  value  problem 

(6.3)  F(x.x')  =  ^  x')  )  ’  x(°)  =  *0.  x'(0)  =  p0- 


If  the  conditions  (Al.2.3)  hold.  then,  for  (xo.po)  in  the  set  EA  deined  by  (4.15).  Theorem 
4.1  ensures  that  (6.3)  has  a  unique  (local)  solution.  Moreover,  on  EA  the  derivative 
DPF  of  the  reduced  mapping  (4.10)  is  non-singular. 

The  nonlinear  system,  to  be  solved  at  each  step  of  the  process,  can  be  written  in 
the  form 

G(x)  = 

where  w  incorporates  all  information  at  the  earlier  computed  points.  All  basic  forms 
of  Newton's  method  are  locally  convergent  if  the  derivative  DG  is  non-singular  at  the 
desired  solution.  Hence  we  are  interested  in  the  non-singularity  of  the  matrix 


{  DFi(x) 

\  DpF2{x,p)  -  hDxF2(x,p) 


DPF{x,p)  +  h 


0  ) 

DxF2(x.p)  ) 


which,  by  definition  of  EA ,  is  clearly  guaranteed  for  all  (x.  p)  in  some  open  neighborhood 
of  any  point  on  the  exact  solution  in  EA  and  for  all  sufficiently  small  h. 

Under  these  conditons  it  can  be  shown  that  when  an  m-step  BDF  method  is  used 
for  the  computational  solution  of  (6.3),  together  with  a  fixed  and  sufficiently  small 
stepsize  h,  then  the  convergence  of  the  approximate  points  to  the  exact  solution  is  of 
order  0(hm)  provided  that  all  initial  points  axe  correct  to  order  0{hm)  and  stopping 
criteria  of  order  0(hm*x)  axe  applied  in  the  Newton  process  at  each  step.  A  proof  of 
this  result  for  the  general  system  (6.1)  may  be  found  in  [3;  where  also  its  extension  to 
the  case  of  variable  steps  is  discussed.  These  results  about  BDF  methods  for  index-one 
systems  form  the  theoretical  basis  for  several  highly  successful  numerical  DAE-solvers 
notably  the  widely  used  codes  DASSL  [30]  and  LSODI  25]. 
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Besides  BDF-methods  also  other  multistep  have  been  considered  in  the  DAE  liter¬ 
ature.  In  particular,  an  extensive  analysis  of  general  linear  multistep  methods  for  the 
index-one  case  is  given  in  22’. 

The  situation  changes  considerably  when  the  DAE  has  index  higher  than  one  and 
two  basic  difficulties  arise.  The  first  derives  from  the  fact  that  for  any  particular 
multistep  method  2  there  exist  DAEs  with  index  exceeding  one  for  which  the  method 
is  unstable.  21.  The  second  difficulty  is  the  appearance  of  a  transient  deterioration  of 
the  discretization  errors  following  any  change  of  the  step-size  in  the  method.  This  type 
of  "boundary  layer”  was  observed  by  several  authors,  see.  e.g..  29.43  . 

More  specifically,  in  43'  it  was  proved  that  when  an  m-step.  fixed- stepsize  BDF 
method  is  applied  to  a  linear  DAE  (4.1)  with  index  k  >  1.  then  the  process  converges 
with  order  0(hm)  after  («  —  l)m  —  1  steps.  Moreover,  in  19  it  was  shown  that  when 
variable  stepsizes  axe  used  and  the  ratios  of  adjacent  steps  remain  bounded  then  the 
global  error  has  order  Oih^^)  where  n  =  min (m.m  —  k  —  2).  Hence,  for  instance, 
for  an  index-three  system  the  use  of  the  implicit  Euler  method  with  variable  stepsizes 
may  lead  to  errors  of  order  0(1).  However,  note  that  for  index-two  problems  we  have 
H  —  m  and.  in  fact,  it  turns  out.  2.4.27],  that  for  semi-implicit  systems  (6.3)  of  index 
two  after  m  —  1-steps  the  m-step  BDF  method  with  fixed  steps  is  globally  convergent 
of  order  0(hm )  provided  again  that  the  initial  points  axe  correct  to  order  0(hm)  and 
the  stopping-criteria  of  the  iterative  process  at  each  step  have  order  0(/im*'1).  But. 
nevertheless,  these  iteretive  methods  may  well  converge  very  poorly  in  the  beginning 
steps.  The  vaxiable-stepsize  case  of  this  result  is  discussed  in  20  . 

Besides  multistep  methods  various  one-step  methods  have  also  been  considered  for 
the  computational  solution  of  DAEs.  In  particular,  there  exists  a  large  literature  on  the 
use  of  implicit  Runge-Kutta  (IRK)  methods.  Any  such  method  can  be  characterized 
by  its  Butcher  tableau 


Cl 

all 

&i2 

•  •  ^lm 

C2 

(221 

022 

•  •  ^2m 

C-rn 

;  ^ml 

Gm2 

-  •  Q-mm 

1  ^ 

..  bm 

see  e.g.  [7].  When  applied  to  the  DAE  (6.1)  the  basic  algorithm  assumes  the  form 

m 

(i)  solve  4-  Cih,xk.x  +  h^di^Yj)  =  0,  z  =  1,2, - m 

j=i 

for  Vj , . . . ,  K,,  £  Rn 

2  In  fact  this  holds  also  for  Runge-Kutta  methods. 
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(n)  set  xk  =  -  h  V 

j=i 


Implicit  Runge-Kutta  methods  are  useful  for  generating  accurate  initial  data  for 
higher  order  multistep  methods;  they  are  also  advantageous  for  problems  with  multiple 
discontinuities.  In  general,  the  nonlinear  system  arising  at  each  step  has  dimension  mn 
and  may  be  very  costly  to  solve  unless  .4  has  special  properties.  Thus  the  complexity 
of  theses  processes  depends  strongly  on  the  form  of  the  coefficient  matrix  .4  =  (au  i. 
Some  important  special  cases  include  the  DIRK-methods  for  which  .4  is  block-lower- 
triangular  with  equal  diagonal  blocks  as  well  as  the  SIRK-methods  where  .4  has  one 
real  eigenvalue. 


In  the  numerical  integration  of  stiff  ODEs  it  has  become  well-known  that  the  com¬ 
puted  solution  often  exhibits  a  disappointingly  low  accuracy  when  compared  with  the 
order  of  consistency  of  the  method.  For  Runge-Kutta  methods  applied  to  a  class  of  stiff 
linear  ODEs  this  was  first  observed  in  -34'  where  it  was  noted  that  for  stiff  problems  the 
order  of  consistency  should  not  be  based  on  the  classical  Lipschitz  condition.  Instead 
in  15;  and  several  subsequent  papers  (see  also  the  monograph  12])  one-sided  Lipschitz 
conditions  were  used  to  introduce  the  concepts  of  B-consistencv  and  B-convergence 
which  provide  order  results  that  correspond  more  closely  to  the  observed  behavior  for 
stiff  ODEs. 

For  an  IRK  method  let  the  stage  order  be  the  largest  integer  r  >  1  such  that  the 
conditions 


j=i 


i  =  1, 2, ....  m, 


are  valid  for  k  =  1, 2, . . . ,  r.  Moreover,  define  the  quadrature  order  as  the  largest  integer 
q  >  1  for  which  the  conditions 


.7  =  1 


k- 1 


hold  for  k  =  1,2 Then  p  =  min(r, q)  is  called  the  internal  stage  order  and  for 
q  >  p  the  classical  nonstiff  ODE  order  p  satisfies  q  >  p  >  p  1.  There  are  examples 
of  stiff  ODEs  and  IRK-methods  where  p  >  p  and  the  observed  order  of  convergence 
equals  the  internal  stage  order  p  (see  e.g.  [12]). 

This  behavior  is  mirrorred  in  the  application  of  IRK  methods  to  DAEs.  In  fact, 
DAEs  have  a  close  relationship  with  stiff  ODEs,  as  is  suggested,  for  instance,  by  the 
fact  that  singulaxly  perturbed  systems 


x'  =  — l(x,  y),  ey'  = /2(x,y) 


with  small  e  >  0  become  a  DAE  for  e  =  0. 
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Consider  again  (6.3)  subject  to  the  conditions  (Al.2.3)  and  with  ( x0,p0 )  €  EA.  For 
the  approximation  of  the  solution  in  EA  we  consider  an  IRK  method  with  a  non-singular 
coefficient  matrix  A  such  that 

1  —  bT A~le,  <1.  e  =  (1.1 . 1)T  €  Rm. 

which  means  that  the  method  is  A-stable  (see  e.g.  7  ).  If  q  >  p  >  1  then  for  a  constant 
(sufficiently  small)  step-size  h  >  0  the  global  error  is  at  least  of  order  0[h^~l )  provided 
any  error  in  the  initial  point  and  the  termination  criterior.  for  Newton's  method  are  of 
the  same  order. 

A  proof  of  this  result  is  given  in  3  where  also  exa  "S  of  problems  are  found 
for  which  the  achieved  orders  are  higher  than  the  stated  -r  bound.  In  essence,  the 
results  also  carry  over  to  semi-implicit  index-two  probler.  see  again  3;  and  also  5’). 

Various  other  one-step  methods  have  been  used  for  me  computational  solution  of 
DAEs.  This  includes,  for  example,  the  Runge-Kutta-Rosenbrock  methods  considered 
in  ’42'  and  the  extrapolation  methods  studied  by  Deuflhard  et  al  (see  e.g.  13'  or  14  ). 
We  shall  not  enter  here  into  any  further  detail. 

6.2  Local  Parametrizations 

In  this  section  we  turn  to  a  local  parametrization  approach  suggested  by  the  existence 
results  of  the  earlier  sections.  It  was  introduced  in  (41]  and  considered  further  in  32.31  . 
and  is  related  to  the  generalized  coordinate  partitioning  technique  used  in  the  numerical 
solution  of  Euler- Lagrange  equations  by  E.  Haug  et  al  (see  (28.48] ) 

Consider  first  the  system  (4.2)  subject  to  the  conditions  ( B1 ,2,3 ) .  and.  more  specifi¬ 
cally.  suppose  that  we  are  in  the  setting  of  Theorem  5.3.  Then,  for  given  (io-  Po-  ~o )  £  A 
we  wish  to  compute  the  (^solution  x  :  J  —  Sx ,  r  :  J  —  S.,  of  the  DAE  (5.17) 
specified  by  b  =  Fi(x0)  that  satisfies  x'(J )  C  Sp  as  well  as  the  initial  conditions 
x(0)  =  io,  i'(0)  =  p0,  z(0)  =  Zq.  For  any  x  6  5X  the  unique  solution  p.z  of  the 
equations 

DF1(x)p  =  0,  F2{x,p,z)  =  0,  (x,p,  z)  6  K 

is  provided  by  the  values  p  =  77(1)  and  z  =  £( x )  of  the  mappings  of  Lemma  5.1.  As  we 
saw,  once  a  procedure  is  available  for  computing  77(1)  and  ((x)  for  any  needed  1  6  Sx , 
the  problem  of  solving  (5.17)  in  Sx  reduces  to  that  of  solving  the  explicit  ODE 

(6.4)  1  =  77(1),  x  G  5*. 

Since  the  desired  solution  1  :  J  — 1 ►  Sx  of  (6.4)  through  Zo  has  to  remain  on  the 
constraint  manifold  Mb  through  z0  it  is  natural  to  work  with  a  local  coordinate  system 
on  Mb .  For  this  we  use  a  simple  class  of  such  coordinate  systems  applied  earlier  in 
other  differential-geometric  numerical  methods  (see  e.g.  [39,40]).  Let  xc  €  Mb  be  any 
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point  on  the  n  -  5-dimensional  constraint  manifold  Mb  through  x0  and  consider  any 
linear  subspace  I  I  R",  with  dimT  =  n  —  r,  such  that 

(6.5)  T~  ~kevDF(xc)  =  {0}. 

If  .4  6  Z(Rn_a.  Rn)  is  any  matrix  with  orthonormal  columns  which  span  T.  then  (6.5 i 
is  equivalent  with  the  assumption  that 

(6.6)  ^  ^Tl{Xc)  )  €  Isom(Rn). 

Hence  there  exist  an  open  neighborhood  Iq  of  the  origin  in  Rn-a  such  that  for  all  u  E  Iq 
the  system 

F^x)  \  -  (  0 
At{x  -  xe)  )  \  u 

has  a  unique  solution  x  =  'P(u)  €  Rn.  It  is  readily  seen  that  the  resulting  mapping  'P 
is  a  C1  diffeomorphism  from  V0  onto  the  relatively  open  neighborhood  'PVq  C  Mb  of  xc. 
This  is  the  desired  local  coordinate  mapping.  Note  that  with  *j(u)  =  ^(u)  —  xc  -  Au 
we  have  AT^>{u)  —  0,  and  thus 

#  :  4o  — *  Rn.  ^{u)  =  xc  —  .4u  —  uj(u)  6  Mb,  u  €  V0. 

shows  that  the  point  x  =  'P(u)  on  Mb  is  obtained  by  adding  to  xc  the  vector  *4u  €  T 
and  the  orthogonal  correction  u j(u)  E  Tx.  The  local  coordinate  system  at  the  point  xc 
of  Mb  is  completely  determined  by  the  matrix  .4  and  hence  we  shall  also  speak  of  the 
local  coordinate  system  induced  by  that  matrix. 

In  practice,  it  is  often  useful  to  work  with  local  coordinate  mappings  defined  by  the 
tangent  space  T  —  ker DF(xc)  at  xc  E  Mb  (see  e.g.  ^40] ).  Suppose  that  we  compute  the 
0  R-  factorizat  ion 

r>Fi(i,)r  =  (<?!.<?,)(  * ) 

where  the  matrices  Q\  E  L(R*,Rn)  and  Q 2  E  X(Rn_’.Rn)  have  orthonormal  columns 
and  R  E  L(R*,R*)  is  upper-triangular  and  nonsingulax.  Then  Q 2  can  be  used  as  the 
basis  matrix  A  of  T  while  the  columns  of  Qi  span  T,'L. 

A  second,  practically  useful  choice  of  a  local-coordinate  space  T  and  its  basis  matrix 
A  consists  in  determining  a  permutation  e^1 ,  e^2, . . . ,  e’"  of  the  standard  basis  of  Rn  such 
that  the  matrix 

A  =  (ei*+l,eJ*+2,...,«r’n) 

satisfies  (6.6).  This  choice  partitions  the  components  of  the  vector  x  =  (xj,  x2, . . . ,  xn) 
into  a  vector  (ij4+1,  xJf+2, . . .  ,Xjn )  of  independent  coordinates  and  the  complementary 
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vector  (a^.x^ . xu)  of  dependent  coordinates.  This  is  the  choice  underlying  the 

mentioned  generalized  coordinate  partitioning  approach  of  E.Haue  et  al  (loc.  cit.). 

Once  a  local  coordinate  system  has  been  chosen  then  it  can  be  shown  (see  41  )  that 
the  ODE  (6.4)  has  the  local  representation 

(6.7)  u  =  ATi?(¥(iO),  u  6  V0  C  Rn". 

This  is  an  (n  - 5)-dimensional  explicit  ODE  without  constraints  to  which  any  standard 
ODE  solver  can  be  applied  as  long  as  the  computed  points  remain  in  Vo- 

This  local  coordinate  approach  can  also  be  carried  over  to  the  second  order  DAEs 
(4.3).  As  before,  suppose  analogously  that  we  are  in  the  setting  of  Theorem  5.3.  Let 
(x0.  y0.  p0.  zq)  —  A  be  a  given  point  for  which  DFi{xo)yQ  =  0  and  suppose  that  we  wish 
to  compute  the  (^-solution  x  :  J  —  Sx,  z  :  J  —  S.  of  (5.17)  specified  by  b  =  FAx 0) 
that  satisfies  x'(J )  C  Sy,  x"(J )  C  Sq  and  the  initial  conditions  x(0)  =  xo,  x'fO)  =  y0- 
x  '(0)  =  q0,  z(0)  =  Zq.  For  any  {x.y)  £  S0  =  Sx  x  Sy  the  unique  solution  q.  z  of 

DFi(x)q-i-  D2FAx)(y,y)  =  0,  F2{x,y,q,z)  =  0,  (x,y,q,z)  £  K 

is  given  by  the  values  q  =  p(x.y)  and  z  =  £(x.t/)  of  the  mappings  of  Lemma  5.3. 
Thus,  when  a  method  for  evaluating  rj(x.y)  and  ( (x.y)  is  available,  then  the  problem 
of  solving  (5.17)  in  S0  is  reduced  to  that  of  solving  the  explicit  first  order  system 

(6-8)  x=y,  y'  =  r)(x,y). 

As  we  know  the  desired  solution  (x,y)  :  J  —  S0  of  (6.10)  through  ( x0,y0 )  remains 
on  the  tangent  bundle  TMb  of  the  constraint  manifold  Mb  through  x0-  Thus  we  have 
to  work  here  with  a  local  coordinate  system  on  TMb ■  As  before,  let  the  matrix  .4  £ 
L(RT’-’,  Fr*)  induce  a  local  coordinate  system  at  the  point  xc  €  Mb  of  Mb .  Then 

(6.9)  0  :  V0  x  Rn~*  — >  TMb ,  Q(u,v )  =  (^(u),  D$(u)u),  u£  V'0,  v  6  Rn_* 

defines  a  local  coordinate  system  on  TMb ■  Ry  restricting  V'o,  if  necessary,  we  can  choose 
some  neighborhood  U0  of  the  origin  of  V0  x  R""*  which  ©  maps  into  E0. 

In  this  local  coordinate  system  the  differential  equations  (6.8)  assume  the  local  form 

(6-10)  u'  =  t-,  v' =  Att,(*(u),D*(u)v),  (u,v)£U0 

and  hence,  once  again  a  standard  ODE  solver  can  be  applied  as  long  as  the  computed 
points  remain  in  Uq 
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6.3  Euler-Lagrange  Equations 

In  this  section  we  sketch  briefly  how  a  multistp  method  might  be  applied  when  the 
local  parametrization  approach  is  used  in  the  numerical  solution  of  the  Euler-Lagrange 
equations  (2.3),  (2.4);  that  is  the  constrained  equations  of  motion 


(6.11)  $(x.t)  =  0 

M(x.t)x"  —  Dx$(x.  t)Tz  =  Q(x.x'.t). 

For  further  detail  we  refer  to  )24.33_. 

In  the  case  of  (6.11)  it  is  easily  seen  that  (C3)  is  equivalent  with  the  assumption 

(6.12)  rankZ)1.<i>(2\  t)  =  s,  yT M{x.t)y  >  0.  Vy  £  kerDx$(x.t) 
which  is  equivalent  with  the  non-singularity  of  the  matrix  of  the  linear  system 


(6.13) 


M(x.t)  Dx$(x.t)T 
Dx$(x.t)  0 


_  l  Q(x.y.t) 
9(x.y.t) 


Thus  under  the  condition  (6.12)  the  general  existence  theory  applies  to  (6.11),  (see  e.g. 

[41])- 

For  given  (x.y.t).  set 


g{x,y.t)  =  ~(D2xx$(x,t)(y.y)  -  D2xt$(x.t)y  -  Dl${x.t) 


in  (6.13)  and  let  q  =  T)(x,y.t ),  z  =  £(x,y.£)  be  the  unique  solution  of  that  system. 
Then  the  problem  of  solving  (6.11)  is  reduced  to  solving  the  explicit  first  order  system 

(6.14)  x'  =  y,  y'  =  r/(x,  y.  t). 


As  in  the  previous  section  let  A  6  L( Rn_*.  Rn)  induce  a  local  coordinate  system  at  the 
current  point  of  Mb  and  introduce  the  corresponding  local  coordinate  system  (6.9)  on 
the  tangent  bundle  TM Then  it  follows  readily  that  the  local  represention  of  (6.14) 
is  given  by 

(6.15)  u  =  v,  v'  =  .4T77('I'(u,  t),  JDu'Ptu,  t)v  D^iu,  t),t). 

Suppose  that  for  its  solution  we  use  a  (consistent)  explicit  multistep  method  of  the 
form 

JTi  m  in  tti 

uk  =  5Z  ajuk-j  +  h  Y,  03vk-r  vk  =  Y,  a:vk-j  ~hY,  3Jvk-j 
j=i  j= i  :=i  j=t 

with  constant  step  h  >  0.  For  the  computation  it  is  advantageous  not  to  work  with 

the  local  variables  u,  v  but  to  transform  all  formulas  immediately  back  to  the  original 

variables  x,  y. 

If  the  approximations  x*_j,  ,yk-j,yl^j,  Zk-j,  j  —  l,2,...,m  of  the  solution  are 
already  available,  then  the  algorithm  for  computing  x*,  t/^,  y'k*Zk  has  the  form: 
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(i)  Set  tk  =  tk_l  *  h\ 

(ii)  Evaluate 


ak  =  -4r{^  ajXk-j  -hf2  3jy^-j)  a'k  =  ATiit  ~h±  3jyLj}- 

J=1  J=1  J=1  J=1 

(iii)  Solve  the  nonlinear  system 

(  $(*.tfc)  \  /  0  \ 

l  ^  J  (a  J 

and  set  xk  =  x: 

( i v )  Solve  the  linear  system 

/  Dr$(xfe.tfc)  \  /  -Dt$(xfc,tfc)  \ 

UT  /'“U  j 

and  set  =  y; 

(v)  Solve  the  linear  system 

/  M(xk,tk)  Dx${xk,t)T  \  /  w  \  /  Q(xk,yk.tk )  \ 

V  0  /\:  /  i  g{xk.yk,tk)  j 

and  set  y'k  —  w  and  zk  =  z. 

In  stage  (ii)  the  multistep  formula  is  evaluated  in  terms  of  the  original  variables 
x,y.  Then  the  stages  (iii)  and  (iv)  determine  the  local  coordinate  mapping  (6.15)  and 
finally  in  stage  (v)  the  linear  system  (6.13)  is  solved  to  obtain  the  accelerations  y'k  and 
the  algebraic  variable  zk  defining  the  constraint  force. 

In  stage  (iii)  a  chord-Newton  process  can  be  used  involving  the  matrix  obtained  in 
stage  (iv)  of  the  previous  solution-step.  When  the  computed  points  leave  the  domain 
of  validity  of  the  current  local  coordinate  system,  then  the  matrix  A  has  to  be  updated. 
The  need  for  this  can  be  detected  by  monitoring  the  number  of  iteration-steps  of  the 
nonlinear  solver  in  stage  (ii)  or  the  condition  of  the  linear  system  in  stage  (iii). 

When  an  implicit  multistep  method  is  to  be  used  then  all  three  steps  (iii)-(v)  have  to 
be  combined  into  one.  Now  special  attention  has  to  be  given  to  the  inherent  struc  .re 
of  the  resulting  large  nonlinear  system  in  order  to  keep  the  computational  complexity 
at  an  acceptable  level.  For  some  detail  we  refer  to  [31]  where  also  a  numerical  example 
is  given. 
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7  Appendix 

In  this  Appendix  we  collect  some  background  material  used  throughout  the  presenta¬ 
tion.  For  further  details,  especially  on  the  differential  geometric  aspects,  we  refer  to 
standard  text  such  as  44'  or  L. 

As  usual,  a  mapping  F  :  U  —  Rm  on  the  open  set  U  Z  Rn  is  of  class  Cr .  r  >  0. 
on  V  if  F  is  continuous  and  for  r  >  0  all  its  partial  derivatives  up  to  and  including 
order  r  exist  and  are  continuous  on  U.  More  generally,  a  map  F  :  S  —  Rm  on  an 
arbitrary  set  5  C  Rn  is  of  class  C"  if  for  each  x  6  5  there  exists  an  open  set  U  I  Rr 
containing  x  and  a  Cr-mapping  F  :  V  —  R™  that  coincides  with  F  throughout  U  ~  S. 
A  map  F  :  5  C  Rn  —  T  Z  R™  is  a  homeomorphism  between  the  sets  5  and  T  if  F 
is  a  one-to-one  mapping  from  5  onto  T  and  both  F  and  its  inverse  F~l  :  T  —  5  are 
continuous.  Finally,  a  map  F  :  5  C  Rn  —  T  C  Rm  is  a  C-diffeomorphism  if  F  is  a 
homeomorphism  between  5  and  T  and  if  both  F  and  F~l  are  of  class  Cr. 

A  subset  M  C  Rn  is  a  d-  dimensional  Cr-  sub- manifold  of  R"  if  for  each  point  x  €  M 
there  exists  an  open  set  V  C  Rn  containing  x  such  that  the  neighborhood  U  ~  M  of  x  on 
M  is  Cr-difFeomorphic  to  an  open  subset  V  of  Rd.  Any  particular  such  diffeomorphism 
O  :  V  n  M  —  V  is  called  a  chart  and  its  inverse  a  local  coordinate  system  on  U  ~  M . 

By  this  definition  any  open  subset  V  Z  Rn  is  an  n -dimensional  Cx -sub-manifold 
of  R".  The  tangent  space  TXU  of  this  manifold  V  at  any  point  x  Z  l  is  defined  as  the 
tt-dimensional  linear  space  {x}  x  Rn.  and  its  tangent  bundle  TV  is  the  2n-dimensional 
submanifold  U  x  Rn  of  R2n. 

Let  F  :  I’  C  R"  x  Rm,  n  >  m,  be  some  C-mapping,  r  >  1.  on  the  open  set  V  Z  Rn- 
A  point  x  €  V  is  a  regular  point  of  F  if  dim  DF(x) Rn  =  m ;  that  is.  if  the  derivative 
DF{x)  has  full  rank  m.  If  all  points  of  a  set  5  C  Rn  are  regular  points  then  F  is  a 
submersion  on  5.  A  point  b  €  Rm  is  a  regular  value  of  F  if  all  points  of  the  inverse 
image  F~l(b)  =  {x  €  U,  F(x)  =  b }  are  regular;  that  is,  if  F  is  a  submersion  on  F~l{b). 
A  fundamental  result  then  states  that  for  any  regular  value  b  £  Rm  the  inverse  image 
Mb  =  F~l(b )  is  either  empty  or  a  p  —  (n  —  m)-dimensional  Cr-sub- manifold  of  Rn.  The 
tangent  space  TxMb  at  any  point  x  of  this  manifold  Mb  may  be  identified  with  the  set 

TxMb  =  {(x,p)  6  Tx Rn;  DF(x)P  =  0}. 

Clearly,  TxMb  is  a  p-dimensional  linear  subspace  of  the  n-dimensional  linear  space 
!TxRn.  The  tangent  bundle  T Mb  of  Mb  is  the  disjoint  union  of  all  tangent  spaces  Tx Mb 
for  x  €  M^  that  is, 

TMb  =  {(x,p)  €  TU;  F(x)  =  6,  DF{x)p  =  0}, 

and  TMb  is  a  2d-dimensional  C^-sub-manifold  of  TRn.  Evidently,  then  the  tangent 
bundle  of  TMb  is  the  4d-dimensional  Cr_2-sub-manifold 

T2Mb  =  {((x}y)(p,q))eT2U;  F(x)  =  b,  DF(x)p  =  0,  DF(x)q  +  D2F(x)(y.p)  =  0}. 
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of  R4n. 

A  Cs-vectorfield  on  some  open  subset  l  C  Rn  is  a  C’-mapping  on  L  such  that 

(7.1)  7r  :  U  —  TV:  tt(x)  =  (1,77(1)),  vx  €  V. 

An  integral  curve  of  ir  through  a  point  x0  £  V  is  any  C'-path  1  :  J  — *  V,  defined  on 
an  open  interval  J  Z  R1  containing  the  origin,  for  which  x(0)  =  x0  and  (if  1 1.  x'l  1 1 1  = 
7r(i(f))  for  t  £  J\  that  is.  which  solves  the  initial  value  problem 

x'  =  77(1),  1  6  i\  x(0)  =  x0. 

For  a  vectorfield  (7.1)  of  class  C‘ ,  s  >  1,  on  the  (non-empty)  open  subset  V  Z  Rn. 
the  following  results  hold: 

(i)  There  exists  a  C*-integral  curve  x  :  J  —+  U  o{  n  through  each  x  £  V  defined  on 
an  open  interval  J.  Moreover,  any  two  such  curves  axe  equal  on  the  intersection 
of  their  domains. 

(ii)  The  union  of  the  domains  of  all  integral  curves  of  n  through  a  point  x  £  V 
is  an  open,  possibly  unbounded  interval  J’.  There  exists  a  C’-integral  curve 
1*  :  J*  — *  U.  of  7r  through  1  and  J’  is  the  largest  interval  on  which  such  an 
intergal  curve  exists. 

(iii)  The  set  D( tt)  =  {(f,x)  E  R1  x  U\  t  £  J*}is  open  in  R1  x  U  and  contains  {0}  v  V. 
Moreover,  the  global  flow  £  :  D(n)  —  U,  £(£.  1)  =  1*,  t  £  J ‘  of  7 r  is  of  class  Ca 
on  D( 7r). 

Consider  now  a  second  order  ODE  x"  —  77(1,1'),  x  £  U  where  77  is  of  class  C* 
on  some  open  set  E  C  R2n.  When  this  problem  is  written  in  the  first  order  form 
x'  =  y,  y'  =  77(1,  y),  (x,y)  €  E,  then  we  encounter  a  vector-field  of  the  form 

(7.2)  ir  :  E  CTU  -7  T2U;  ir  (z,y)  =  ((1,  y),  (y,77(i,y))),  V(i,y)  6  E. 

Note  that  the  second  and  third  component  of  the  image  vector  are  identical;  in  other 
words,  (7.2)  represents  a  sub-class  of  all  tangential  vector  fields  on  TV,  namely,  the 
vector  fields  that  are  consistent  with  the  second  order  ODE. 

An  integral  curve  of  (7.2)  through  a  point  (i0,J/o)  €  E  is  now  a  C'-path  x  :  J  —  Rn. 
defined  on  some  open  interval  J  C  R1  containing  the  origin,  for  which  (x(t),x'(t))  €  E 
and  ((x(t),x'(t)),  (x'(t),  x"(t)))  =  7r((x(t),i'(t))  for  all  t  £  J;  that  is,  which  is  a  solution 
of  the  original  second  order  ODE. 
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